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ABSTRACT 

We investigate the vertical structure of neutrino dominated accretion disks by self- 
consistently considering the detailed microphysics, such as the neutrino transport, ver- 
tical hydrostatic equilibrium, the conservation of lepton number, as well as the balance 
between neutrino cooling, advection cooling and viscosity heating. After obtaining the 
emitting spectra of neutrinos and antineutrinos by solving the one dimensional Boltz- 
mann equation of neutrino and antineutrino transport in the disk, we calculate the 
neutrino/antineutrino luminosity and their annihilation luminosity. We find that the 
total neutrino and antineutrino luminosity is about 10 54 ergs/s and their annihilation 
luminosity is about 5 x 10 51 ergs/s with an extreme accretion rate 10M sun /s and an 
alpha viscosity a = 0.1. In addition, we find that the annihilation luminosity is sensitive 
to the accretion rate and will not exceed 10 50 ergs/s which is not sufficient to power 
the most fireball of GRBs, if the accretion rate is lower than lM sun /s. Therefore, the 
effects of the spin of black hole or/and the magnetic field in the accretion flow might 
be introduced to power the central engine of GRBs. 

Subject headings: accretion, accretion disks, black hole physics, gamma ray bursts, 
neutrino 



1. Introduction 

Gamma ray burst(GRBs) are e xtremely high energy rel e asing phenomena in the universe and 



are usually divided into two classes (jKouveliotou et al.lll993l ; IZhang fe Meszarosl 12004 : iPiran 



2004 



Nakarl 120071 ) : short GRBs (T90 < 2s) and long GRBs (T90 > 2s). Numerous models have been 
proposed to explore the central engines of GRBs, and one of the mostly discussed model is the 
neutrino dominated accretion flows (NDAFs) with a hyper accreting stellar massive black hole 
with accretion rate 0.1 ~ 10M sun /s. Due to the high density and high temperature in the inner 
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part of NDAFs, the optical depth of photons is very large and photons are completely trapped, 
then neutrinos and antineutrinos become the most promising candidates that carry away thermal 
energy and cool the disk. The annihilation of n eutrino pairs above the disk is believed to be 
the energy source of GRBs. iNaravan et al.l (|1992l ) first proposed that neutrino pairs annihilation 
into electron pairs d uring the merger of compact objects binaries may power GRBs. After that, 



Popham et al.l dl99J) investigated NDAFs under the assumption that the disk is transparent to 
neutrinos, but they also pointed out that the assumption fails when the accretion rate is higher 
than lM sun /s . |Pi Matteo et al.l (|2002i ) improved the model by using a simplified neutrino transport 



model which was bel i eved to bridge the neutrino optically thin limit and optically thick limit. 



Chen &: Beloborodovl (|2007l ) improved the model further by dealing with the neutrino emission 



and chemical composition in the optically thin regime, optically thick regime and intermediate 
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there are still some uncertainties including: 1, The distribution of electron fraction. In many 
previous works, which is assumed to be a constant value, for example, 0.5, throughout the disk. 
While th e electron fraction has a g r eat effect on the emission of neutrin os and antineutrinos, as 
shown by iKohri fc Mineshigd ([2002) ; iKohri et al.l (|2005l ) ; iLiu et al.1 (120071 ) , so need more cautious 
disposal. 2, Neutrino transport and neutrino spectra. The most commo n ly use d approximation was 
the simplified neutrino transport model introduced in lDi Matteo et al.l ((2002). In their model, the 
di fference between the neutrino and antineutrino transport in the disk was neglected, but as shown 
in lPan &: Yuan! (|2012j ). the precise spectra of neutrino and antineutrino sensitively determines the 
annihilation luminosity of neutrino pairs. 3, The annihilation of neutrino pairs . The most common 
method for calculating the annihilation luminosity was originally introduced by Ruffert et al. (jl997l ) 
to calculate the annihilation luminosity during the merger of neutron star binaries. This method 
was applied in the calculation of the annihilation luminosity above NDAFs under the assumptio n 
that the emission of neutrinos and antineutrinos are isotropic and symmetric (jPopham et al.lll999l ). 



It is evident that the most strict approach to determine the neutrino/antineutrino luminosity 
and their annihilation luminosity above NDAFs is to build the two-dimensional disk model in which 
neutrino transport, vertical structure, chemical evolution, thermal evolution and the distribution 
of electron fraction, mass density, and temperature are self consistently considered. 



Rossi et al.l (|2007l ) first investigated the vertical structure of NDAFs by using the Eddington 
approximation to deal with the neutrino transport in the vertical direction, neglecting the con- 
tribution of advection term to the disk cool ing, and deal i ng wi th ne utrino emission and chemica l 
composition following the similar method of I Janiuk et al.l (|2007l ) and lChen &: Beloborodovl (|2007l ). 



Recently, ILiu et al.1 (120101 ) also investigated the vertical structure of NDAFs with many sim- 
plifications on neutrino transport, equation of state and the annihilation efficiency. The first one is 
that the neutrino emission is directly integrated to calculate the neutrino luminosit y by neglecting 



the absorption of neutrinos, which is viable in the neutrino optically thin limit (jPopham et al 
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1999 ); the second one, that a simplified equation of state p = Kp 4 ^ 3 is used in the vertical direc- 
tion, which is viable when relativistic degenerate electrons dominate the pressure of the disk; the 
third one that a toy annihilation efficiency of neutrino pairs r\ = L uD /L u oc V^j* is applied, where 
L V v,L v is the annihilation luminosity and neutrino luminosity before annihilation respectively, and 
Vann is the so called annihilation volume. With all the above simplifications, the annihilation effi- 
ciency r\ = L V yjL v can even reach 100%! Obviously, the unrealistic result is caused by too many 
unrealistic assumptions. 

In this work, we investigate the vertical structure of NDAFs, neutrino/antineutrino luminosity 
and their annihilation luminosity by self-consistently considering the neutrino/antineutrino trans- 
port, the vertical hydrostatic equilibrium, the precise equation of state, the chemical equilibrium 
and the thermal balance between neutrino cooling, advection cooling and the viscosity heating 
under the self similar assumption of the distribution of mass density and i nternal energy density in 
the radial direction (jShakura k, Sunvaevlll973l ; iNaravan &: Yilll994l . Il995l ). Especially, we strictly 
solve the Bolzmann equat ion to deal with the neutrino transport, instead of taki ng the assumptio n 
of the gray body spectra (IJaniuk et al.1 120071 ) or the Eddington approximation (IRossi et al.ll2007l ). 
Correspondingly, we can precisely obtain the energy spectra of neutrino pairs. Combining the con- 
servation of the number of lepton, the distribution of chemical composi tions are self-consistently 
and accurately determined (jChen & Beloborodovl 120071 : IRossi et al.l 120071 ). 



This paper is organized as follows. In §11, we introduce the basic equations in our calculation, 
including the Boltzmann equation of neutrino/antineutrino transport, angular momentum equa- 
tion, hydrostatic equilibrium equation, equation of state, lepton number conservation equation and 
thermal evolution equation. In §111, we briefly introduce our numerical methods to find the steady 
solution of the structure of the disk. In §IV, we list our numerical results of the disk structure, 
neutrino/antineutrino luminosity and their annihilation luminosity. Conclusions and discussions 
are summarized in §V. 



2. Basic Equations 



We assume a steady accretion disk with accretion rate M = 10, 1, 0.1M sun /s around a cen- 
tral black hole with mas s M = 3.3M sun and adopt the standard a viscosity prescription of 
Shakura k, Sunyaevl ()1973l ) with a = 0.1 for the viscous stress of the disk. We discuss the struc- 
ture of the disk and neutrino transport in cylindrical coordinate (r, z, (j)) and we assume the inner 
boundary of the disk to be r- m = 6M and outer boundary to be r ou t = 100M. 



2.1. Boltzmann equation 



We solve the one dimensional Boltzmann equation of neutrino and antineutrino transport in 
the vertical direction of the disk and obtain the energy dependent and direction dependent neutrino 
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spectrum. We define //) and ^x) to be the distribution function for up moving neutri- 

nos /antineutrinos and down moving ones respectively, where z is the vertical coordinate of the disk, 
p is the energy of neutrinos/antineutrinos, and p = cos(#) for up moving neutrinos/antineutrinos 
and p = — cos(#) for down moving ones, where 6 is the angle of neutrinos/antineutrinos moving 
direction to the vertical direction of t he disk. For the up- moving neutrinos/antineutrinos, their dis - 



tribution function is determined by (ISawyerll2003l : iBurrows et al. 

df+(z,p,p) 



2006 



Schinder fc Shapirolll98^ : 



/'- 



dz 



\ a [f*(T(z),p eq ,p) ~ f+(z,p,p)]+\ s 



-f+{z,P,p) + 



1 



dfJ>f-(z,p, p) + f+(z,p,n) 
(1) 



and for the down-moving ones, their distribution function is determined by 
df-(z,p,fi) ! /J 



P 



dz 



-A„ [f eq (T(z),p eq ,p) - f-(z,p,p,)]-X s 



/_(z,p,/i) + i J dfif-(z,p,fi) + f+(z,p,fi) 

(2) 



Where A a is the absorption coefficient and A s is the scattering coefficient of neutrinos/antineutrinos, 
and here / eq = l/(exp ((p — fi eq )/kT) + 1) for neutrinos, / eq = l/(exp ((p + fx eq )/kT) + 1) for anti- 
neutrinos, where fi eq = p c + /i p — /i n , and fj, e , yu p , fx n is the chemical potential of electron, proton 
and neutron, respectively. 

Because Urea process v e + n o e~ + p and v e + p o e + + n dominate the creation and 
the absorption of neutrinos and antineutrinos, and the neutrino/antineutrino scattering by neu- 



trons and proto ns (V fi , v e ) + p 



Ve) + p and u P ) 



n 



{ v e , + n dominates the scat- 



tering opacity (jPopham et al.l [l999j; iJaniuk et al.l 120071 ; iLiu et al.l 120071 ) , so we include no other 
neutrino/antineutrino processes. Thus, for simplicity, in this paper we use v and f e , v and v e in- 
terchangeably and it will not cause any confusion. As for the explicit expression of t he absorption 
coeffic ient A a and scattering coefficient A s of neutrinos/antineutrinos, please refer to |Pan Yuan 
|2Qll). 



Considering that the disk is symmetric about the equator plane, the boundary conditions 
for the distribution function f(z,p,fi) of neutrinos/antineutrinos can be written as f + (0,p,fi) = 
/L(0,p, yu) and f-(H,p,fi) = 0, where H is the upper boundary of the disk. 



2.2. Angular momentum equation and Vertical hydrostatic equilibrium equation 

Adopting a prescription, the tangential stress and angular momentum equation can be written 
as (jShakura &: Sunyaevlll973l ). 



Wrcf, = apc s , 



(3) 



and 



d(Qr 2 
1 dt 



pv, 



d(flr 2 ) _ 1 d{w r(j) r 2 ) 



dr 



dr 



(4) 
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where w r( p is the viscous stress, p is the mass density of the disk, c s is the acoustic speed, v r is the 
radial drift velocity and fi is Kepler angular velocity at radius r. Integrating Eq.(4) over radius r 
and height z, we obtain 

f H f r d(flr 2 ) , , f H f r d(w r ^r 2 ) , , 

2vr / / rpv r ^——^dzdr = 2vr / / v r f ' dzdr, (5) 



HJr in dr ' J-HJr in dr 

and by taking into consideration the steady accretion condition 

M = 2nr / pv r dz = const, (6) 

J-H 

and using torsion condition in the inner boundary ttv^( r m) = 0, the angular momentum equation 
Eq.(5) is transformed to be 



M (Qr 2 - (fir 2 ) in ) = 2vrr 2 / apc 2 dz. 

J-H 



The vertical hydrostatic equilibrium equation is simply as follows, 

pGM z _ dp 
r 2 r dz 



(7) 



(8) 



2.3. Equation of state (EOS) 



In our calculation, the EOS of accreted gas including protons, neutrons, and electron pairs are 
determined by the exact Fermi-Dirac integral, the EOS of neutrinos pairs are determined by the 
numerical integration of their distribution functions [f+-(z,p, p)] u ,u, and the EOS of photons is 
simply according to Eq.(15) and we do not include other kinds of particles in this work, especially, 
helium, which was included in the most of the previous works (we will check the validity of neglecting 
the contribution of helium in §IV): 



p(p, Y C ,T) = p n +p p +p e + p e + + p rad + p u + p p , 

U(p, Y e , T) = U n + Up + lie + u e + + 1i rad + u v + Up, 



(9) 
(10) 



where p and u is the total pressure and total internal energy density respectively, Y e is the electron 
fraction Y e = (n e — n e +)/(n p + n n ), and T is the local temperature of the disk. 



Specifically, the EOS of gas are expressed as (jJaniuk et aJj 120071 ). 



Pi 



2^2 Kc 2 ) 4 5/2 



3vr 2 (he) 



-ft 



1 



^3/2(r?i,ft) + -A^5/ 2 (^,ft) 



^ W^ /2 + ^s/ 2 (»»,A)] 



rriiC 



IT 2 \ he 



Pi /2 [F 1/2 (rii, Pi) +m/2(m, Pi)] 



(11) 

(12) 
(13) 
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where pi, Ui,rii is the pressure, internal energy density and number density of particle i respectively 
(i = n, p, e, e + ), F^ is the Fermi-Dirac integral of order k, rji is the degeneracy parameter of particle 
i (rji = tif /kT, where nf is the chemical potential of particle i not including the rest mass) and 
fii is the relativity parameter of particle i = kT/rriiC 2 ). 

And the EOS of neutrino pairs and radiation are expressed as 

Pi = f , (14) 



vr 2 (kT) 4 
15 [hcf 



■"rad = — 7i s 3 , (15) 



Uu,u = 0- J J P 3 (f+ + f-)u,udpd/i, (16) 

n u ,u = ~f~§ff p2 (f+ + f-)v,vdpdii. (17) 

where p is the energy of neutrinos and antineutrinos and pj is the pressure of particle j (j 
rati, v, v). 



2.4. Lepton number conservation 

It is easy to write down the equation of lepton number conservation of fluid using the Eulerian 
description: 

d(n h Y lep ) dF lep 

+ v • V(n b Yi ep ) + = 0, (18) 

where rib is the number density of baryon, Y\ ep is the fraction of lepton number 

v _ n e - ri e + + re^ - n F 

-Men — j V iy J 

n b 

and Fi e p is the lepton number flux in the vertical direction of the disk ( here F\ ep is contributed by 
neutrinos and antineutrinos F\ ep = F u + Fp) 



2vrc / ; 2 



F » = -^J J P (f+ ~ f-),ufidpdfJi, (2°) 
Fp = ~ [ J P 2 (f+ - f-),^dpdfj,. (21) 

Due to the fact that the time scale of accretion is much longer than that of chemical evolution 
under the condition of the inner part of NDAFs: p « 10 10 ~ 10 11 g/cm 3 , or rib ~ 10~ 5 ~ 10" 4 
fm -3 , and T ~ 5 x 10 10 K, the time sca le for the typical Urea process p + e~ — > n + v e is about 



0.1 ~ 1 ms ( see the Fig. 2 of lYuanl ([2005) ) which is much shorter than the time scale of accretion, 
so it is reasonable to neglect the advection term in the equation of lepton number conservation, 
hence Eq.(18) is simplified to be 

at + dz " u - 1 ] 
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2.5. Thermal evolution equation 



The energy equation of fluid is written as 

du 



dt 



+ v • Vu 



pV • v, 



or in the equivalent form 



du 
dt 



Q+ - q- 



Qadv 7 



(23) 



(24) 



where q a ^ v = pV • v + v • Vu is the advection cooling term and q+ is the alpha-viscosity heating 
rate 



q+ = apc s r— 



dr J ' 



(25) 



and q- is neutrino cooling rate contributed by the energy flux of neutrinos and antineutrinos: 
lu + Qu 

P 3 (f+ - f-)^dpdu. (26) 



2irc d 



h 3 dz 



We assume the velocity distribution to be v r = —acl/rQ, v z = 0, and U = y GM/r 3 , which 
is similar to the self similar radial d istribution of the mass density of the gas pressure dominated 
thin disk (jShakura k, Sunvaevill973l ). we also assume the distribution of mass density and energy 



density to be self similar p 



~ 3 I 2 and 



u ~ r 



-2 



in the radial direction. According to the mass 



conservation equation of steady flows pV • v + v • Vp = 0, we obtain V • v = 3/2(v r /r), so the 
advection term q^ v is simplified to be 



iadv 



7 r 2n 



(27) 



and we will check the validity of the self similar assumption in §IV. 



For simplicity, we adopt the value of acoustic speed at z = when calculating viscosity 
heating rate and radial velocity at any location, i.e. more numerically economic expressions w r ^ = 
ap [c 2 s (z = 0)] for the viscous stress and v r = a\c 2 ( z = 0)1 /rVt for the radial velocity, and c 2 . is the 
square of adiabatic sound speed, which is given by (|Fontl 120031 ) 



dp 



d(p + u) 



1 



ad 



p + p + u 



Op 



dp 



dp 



du 



(28) 



3. Numerical Methods 

3.1. Two stream approximation 

The two-stream approximation is a simplification to the full Boltzmann equation that replaces 
the full direction dependent distribution f + (z,p,p) and f-(z,p,p) by two streams with angle 
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cos(#) = ±l/y3 to the vertical direction. Under the two stream approximation, the full Boltzmann 
equation is simplified to be 

1 df+(z,p) 



V3 dz 
1 df-(z,p) 



\ a [f*(T(z),p cq ,p) - f + (z,p)] + -X s [f-(z,p) - f+(z,p)], 

~ A a [/ 6q (T(z), /i eq , P) ~ f-(z, P )] + \\s\f-{*,P) ~ /+(*,2>)]- 



(29) 

Two stream approximation has been confirmed to be valid by Sawyer (j2003l ): Pan Sz Yuan ( 2012 ) 
and is much easier to deal with compared with the full Boltzmann equation, so it is a prime choice 
for the calculation of the energy density u Vt p, the cooling rate q U) „ and the lepton number flux F V) p 
of neutrinos and antineutrinos: 



F - 



1 2vrc d 
V3 h 3 dz 
1 2vrc 



V3 h? 



2ir f 

^3 / P 3 (f+ + f-)u,udp, 

P 3 (U ~ f-)u,udp, 
P 2 (f+ - f-)u,vdp. 



(31) 
(32) 
(33) 



3.2. Numerical methods 

Now, we start to seek for a steady solution to the disk in which hydrostatic equilibrium, 
thermal balance and chemical balance are satisfied. First, we assume an initial distribution of 
temperature T°(r, z) = 8(r/r; n ) 1//3 MeV and election fraction Y®(r, z) = 0.4 (In fact, our numerical 
calculation has shown that the final convergent solution does not depend on the specific choice of 
the initial conditions which only affects the speed of convergence of numerical calculation. Such 
independence proves the validity of our numerical methods in turn), then we solve the vertical 
hydrostatic Eq.(8) subjected the boundary condition (7), and it is not hard to obtain the zero 
order mass density distribution p°(r, z) and corresponding the thickness of the disk H°(r). Then, 
we solve the two stream approximation Eq.(29) and (30), and obtain the neutrino/antineutrino 
spectra f±^ uj) y With the neutrino/antineutrino spectra, we can solve the chemical evolution Eq.(22) 
and thermal evolution Eq.(24) with the zero order mass density distribution p°(r, z) or the baryon 
number density n®(r,z) fixed , until a steady state dT l /dt = and dY^/dt = 0. With the first 
order distribution T x (r, z) and YJ-(r, z), we solve the vertical hydrostatic Eq.(8) subject to the 
boundary condition (7) once more to get the corresponding first order mass density distribution 
p l {r,z). Iterating this process until the final overall convergent solution of mass density p(r,z), 
electron fraction Y e (r, z) and temperature T(r,z) is obtained. 

After that, we switch to the full Boltzmann equation Eq.(l) and (2) to solve the full en- 
ergy dependent and direction dependent distribution function of neutrinos and antineutrinos: 
[f +i -(z,p, and [f + ^(z,p, p)] t p. With the full energy dependent and direction dependent spec- 
tra of neutrinos and antineutrinos, we can calculate their annihilation luminosity precisely. The 
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annihilation rate of neutrino pairs is expressed as follows (jRuffert et al.lll997l ; iBurrows et al.ll2006l ): 



Q{v e v e 



cr c 



C X +C 2 



pOO pOO p p 

I dp dp'(p + p')(pp'f dn dn 1 'f v jv.{i - ^se) 5 

J0 JAtt Ji-K 

.(34) 



4 (m c c 2 ) 2 (hcf 

poo poo p p 

+ c 3 (m e c 2 ) 2 dp dp'ip + p'^pp') 2 dn dn'f u j Vc (i- cos e) 

Jo Jo J4n JAtt 



where the typical cross section of neutrino interaction is do = 1.705 x 10 -44 cm 2 , the weak interaction 
constants are C\ + C 2 ~ 2.34, C3 ~ 1.06, p and p' is the energy of neutrinos and antineutrinos 
respectively, n and n' is the solid angle of the incident direction of neutrinos and antineutrinos 
respectively, O is the angle between neutrino beams and antineutrino beams (see Fig.l). 



4. Results 

In this section, taking the case M = 10M sun /s, a = 0.1 as an example, we show the numerical 
results of the vertical structure of the accretion flow, its radial structure, the spectral energy distri- 
tion of the neutrino/anti-neutrino from the disk and the final annihilation luminosity of neutrinos 
and anti-neutrinos. 



4.1. The vertical structure of the disk 

Figure 2 shows the distribution of the mass density p, temperature T, electron fraction Y e and 
the ratio of internal energy density to pressure u/p in the vertical direction at two different radius 
10M and 35M, respectively. 

The mass density is about p ~ 10 11 g/cm 3 and the temperature is about T ~ 8 MeV in 
the inner part (r = 10M) of the disk. The precise value of the temperature in the inner part of 
NDAFs is vitally important, which sensitively determines the annihilation luminosity of neutrino 
pairs L V y ~ T 9 according to Eq.(34). 

When the disk is in chemical equilibrium, the electron fraction cannot be describe by only a 
constant value throughout the disk, while it varies a few times from the bottom to the surface of the 
disk. It is the specific distribution of the electron fraction that guarantees the chemical equilibrium, 
which was not included in the most previous works, as a result, the electron fraction has to be an 
artificial assumption there. 

It is noticeable that the vertical distribution of u/p and electron fraction Y e are positively 
correlated (Fig. 2c and 2d). The positive correlation implies that relativistic electrons dominate 
the pressure at the surface of the disk where u/p ~ 3, and non-relativistic protons and neutrons 
dominated the pressure at the bottom of the disk where u/p ~ 3/2. To justify the conclusion, 
we plot the vertical distribution of the ratio pi/p at radius 10M and 35M in the Fig. 3, where 




Fig. 2. — The vertical structure of the disk at different radius: r = 10M (solid lines) and r = 35M 
(dashed lines), where H is the thickness of the disk (see Fig. 6) 
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Pi = Pp +Pn, P2 = Pe + P e + an d P3 = Prad +Pf It is indeed so, baryons (neutrons and protons) 
dominate the pressure at the bottom of disk and leptons (electron pairs) dominated the pressure 
at the surface o f the disk. So it is not reasonable to assume a polytropic equation of state p oc p 4 / 3 
(|Liu et aL||2Q10D which is the EOS of relativistic and strongly degenerate electron gas whose ratio 
of internal energy density to pressure is u/p = 3. If a polytropic equation of state is needed to do 
some approximation and estimation, actually p oc p 5 / 3 is a better choice. 



4.2. The radial structure of the disk 

Two main assumptions are used in the above discussion: we apply the self similar assumption 
of the distribution of mass density p and internal energy density u in the radial direction when 
calculating the advection cooling term and we neglect the contribution of helium to the EOS and 
the contribution of helium disintegration to the cooling term. We now check their validity. 



4-2.1. Self-similar behavior in the radial direction 

Fig. 4 shows the radial distribution of mass density p, internal energy density u, pressure p and 
temperature T on the equator plane z = and their corresponding fitting lines: p ~ r _1 ' 5 ,u,p ~ r~ 2 
and T ~ r ~°- 6 . According to Fig. 4, the self similar assumption of mass density p and internal energy 
density u in the radial direction is rather a self consistent and accurate description. 

Now we give a more physical explanation about the self-similar behavior in the radial direction: 
according to Fig. 2 and Fig. 3, it is evident that non-relativistic protons and neutrons dominate total 
pressure at the bottom of the disk and determine the surface density of the disk. Hence, we use the 
polytropic equation of state p oc p 5 / 3 to estimate vertical structure of the disk. Combining Eq.(7) 
and (8), it is easy to get the scale thickness of the disk H ~ c s /£l, and the mass density p ~ r~ 3 / 2 . 
In order to estimate the radial behavior of temperature and pressure, we must consider the more 
general EOS of non-relativistic gas p = pT, and the thermal balance between viscosity heating and 
neutrino cooling q+H ~ T 4 or pVLH ~ T 4 , so T 4 ~ pc s , consequently T ~ r~ ° -6 and p ~ r~ 2,1 (see 
Fig4). 

In order to gain an insight to the nature of the disk we investigate, we also calculate the 
advection factor 

/adv = / Qadvdz/ / q + dz, (35) 

where (; a dv and g+ are the advection cooling rate and heating rate defined in the §11, and the 
result is shown in Fig. 5a. So it is evident that the disk is indeed a neutrino cooling dominated 
accretion disk as its name suggests: neutrino radiation dominates the cooling process and advection 
is always a minor role. The self similar assumption is only applied in the calculation of the advection 
term, so even if there is some small deviation between the realistic distribution and the self-similar 
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assumption as shown in Fig. 4, it has no much influence on the final results. Thus the self similar 
assumption is rather a reasonable simplification. 



4-2.2. The fraction of He 



The number density of helium nn e is expressed as 

1 



-(i-x B 



(36) 



where the fraction of free nucleons X nnc (protons and neutrons) is given by (jJaniuk et al 
Qian fe Wooslevlll996h 



2007 



X nuc = 295.5p 10 3/4 r i 9 / 8 exp (-0.8209/T n ). 



(37) 



where pio is the mass density in unit of 10 g/cm and T\\ is the temperature in unit of 10 K 



and if l„„ r > 1, then X n 



1. 



We plot the free nucleons fraction of the gas on the equator plane log(X nuc ) calculated from 
Eq.(37) versus radius log (r/M) in Fig. 5b. According to Fig. 5b, it is easy to know that X nuc 3> 1, 
i.e., the assumption X nuc = 1 and nn c = perfectly hold here. So it is sound to neglect the 
contribution of helium to EOS and the contribution of the helium disintegration to disk cooling. 

It should be noted that the radius where Helium dissociation becomes important depends on 
the viscous parameter a lChen fe Beloborodovl (|2007l ) . So in the case of smaller a, the contribution 
of helium may not be negligible. 



4.2.3. Profile of the disk 

In addition, we plot the radial distribution of the surface density a in Fig. 6a and the profile of 
the disk thickness H(r) in Fig. 6b. Here, the thickness H(r) is defined according to the mass density 
contrast p(z = H) = p(z = 0)/100. H(r) serves as the upper boundary of Boltzman equation (see 
§2.1), and it is different from definition of the usual scale thickness of the disk which is mostly 
used in the one-dimensional disk model. According to Fig. 6, the disk we consider is a geometrically 
thick disk with thickness H ~ r, which will contribute to higher annihilation efficiency as we will 
discuss in the next section. 



4.3. Neutrino spectrum 

When the disk is in chemical equilibrium, the flux of lepton number F\ ep = F v + Fp vanishes, 

J P 2 (f+ ~ f-)iv pdpdp = J p 2 (f + - fidpdp. (38) 
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Fig. 5. — (a)Left Panel: The radial variation of the advection factor / a d v ; {h)Right Panel: The 
radial distribution of free nucleons fraction log(X nuc ) calculated from Eq.(37). 




Fig. 6. — (&)Left panel:The radial distribution of surface density a in unit of 10 17 g/cm 2 ; (h)Right 
panel: The profile of the disk thickness H(r) versus radius r. 
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Figure.7 shows the direction averaged neutrino/antineutrino spectrum 



F 



num — 



(39) 



at the surface of the disk at different radius, where / = f + (H,p, fi), T = 7.4 MeV for r = 10M and 
T = 3.5 MeV for r = 35M (see Fig.2b). Then, at the surface of the disk the chemical equilibrium 
condition Eq.(38) is transformed to be 



So, chemical equilibrium only requires that the total amount of neutrino number flux and antineu- 
trino number flux are equal, i.e. integration of F num of neutrinos and neutrinos over energy p are 
equal, while it is interesting that the form -F num of neutrinos and antineutrinos almost coincide with 
each other in fact. 

In addition, it is necessary to explain that there is a cusp on the antineutrino spectrum in the 
lower energy end: it is easy to know that there is a lower energy limit -Emm = m e + + m n — m p for 
antineutrinos from the Urea process V e +p^-e + +n which gives rise to the cusp on the energy 
spectrum, while there is no such an energy limit for neutrinos from the process v e + n <H> p + e, so 
the neutrino energy spectrum is smooth as expected. 



In order to calculate the neutrino/antineutrino luminosity and their corresponding luminosity, 
we have to make a discount on the vertical structure of the disk: we assume the disk to be lying on 
the equator plane, and then modify the resulting annihilation luminosity by taking the thickness of 
the disk into consideration. With the thin disk simplification, all the elements in the Eq.(34) for the 
annihilation rate are available (see Fig.l), so it is easy to numerically do the integration of Eq.(34) 
over the entire surface of the disk and the whole energy span of neutrinos and antineutrinos. 

In addition, the neutrino/antineutrino energy flux F Ut p, luminosity at the surface of the disk 
FuyLp and their annihilation luminosity L^p are expressed as 




(40) 



4.4. Annihilation luminosity 




(41) 



(42) 



(43) 



(44) 



-16- 



where r- m = 6M, r out = 100M in our calculation and H(r) is the thickness of the disk at radius r 
(see Fig. 6b). 

The results are as follows: L v « L v = 5.2 x 10 53 ergs/s, L vV = 1.66 x 10 51 ergs/s, the 
annihilation efficiency rj = L u p/(L U + hp) = 0.16%. While it is not the final result, we must take 
into consideration of the thickness of disk H ~ r (see Fig. 6b) to correct the thin disk simplification. 

Fig. 8 shows the concrete distribution of the annihilation rate Q u p calculated with the thin disk 
simplification in the r, z plane. The distribution of the annihilation rate is nearly isotropic aside 
the directions occupied by the disk and the half open angle of empty funnel along the central axis 
of the disk is about 45 degrees which is determined by the ratio of thickness to radius of the disk 
H/r ~ 1. The major contribution to the total annihilation luminosity is concentrated at the zone 
near the central black hole of the disk. The solid angle that the assumed thin disk lying on the 
equator plane subtends to the zone is about Othin ~ 2-7T, while solid angle that the realistic thick 
disk subtends to that zone is about f2 t hj c k; ~ 27r(l + cos(45°)) ~ 1.70, and so final annihilation 
luminosity should be multiplied by a modification factor (^thick/^thin) 2 ~ 3. 

Hence, the final results are as follows: neutrino/antineutrino luminosity is about L u = Lp = 
5.2 x 10 53 ergs/s, annihilation luminosity is about L v p = 5 x 10 51 ergs/s and annihilation efficiency 
is about r] = 0.48%. 



5. Conclusions and Discussions 



Temperature of the inner part of NDAFs sensitively determines the final luminosity of neutrino 



10M sun/s and 



annihilation as L u p oc T 9 . In the extreme model of this paper (M = 3.3M sun , M 
a = 0.1), it is found that temperature in the inner part of the disk is about T ~ 8 MeV, and 
the corresponding annihilation luminosity is about L v p ~ 5 x 10 51 ergs/s which roughly satisfies 
the energy demand of the most energetic GRBs. If the temperature of disks changes from 8 MeV 
to 4MeV, the annihilation luminosity reduces about 500 times, which is not enough to power 
the fireball of the most energetic GRBs with the isotropic luminosity about 10 52 ergs/s. In the 
standard model of GRBs, the duration of the energy injection is the duration of the GRB prompt 
emission, which could be about 2s for short GRBs and about 100-1000s for long bursts. Assuming an 
extreme constant accretion rate of 10 Msun/s means therefore an enormous total mass consumption 
for powering GRBs, especially for powering long bursts, which might be realistic. Therefore, the 
simple NDAF model which is investigated in this work can not produce sufficient energy to power 
GRBs, the effects of the spin of black hole or/and the magne t ic field in the accretion flow might be 
introduced to power the central e ngine of GRBs (see 
20091 . bold : l,Taniuk fc Yuanlhoiol . for instance). 



Li 



2002; iGhen k, Beloborodovil2007l : iLei et al 



It is very easy to understand that the temperature of the inner disk mainly depends on the 
accretion rate of the flow. Besides the extreme model, we also investigate the models with the 
lower accretion rate, such as M = 0.1M sun /s and lM sun /s, and the results including the resulting 
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Fig. 8. — Contour of annihilation rate log Q v p/ (ergs- cm -s ) in r, z plane, and the blank zones 
are occupied by the disk. 
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temperature T in the inner part of disk (r = 10M), neutrino and antineutrino luminosity L u + Lp, 
annihilation luminosity L uu and the corresponding annihilation efficiency ij are listed in Table 1. 
For all the three different accretion rate 10, l,0.1M sun /s, the thickness H(r) we defined in §4.2.3 is 
roughly equal to radius r. 

It is obvious that, accretion rate sensitively determines the annihilation luminosity of neutrino 
pairs and the annihilation luminosity will not exceed 10 50 ergs/s when the corresponding accretion 
rate is lower than lM sun /s. So NDAFs with accretion rate lower than lM sun /s are unlikely to serve 
as central engines of GRBs. 

Neutrino and antineutrino spectra is the second major factor determining the annihilation 
luminosity. The resulting neutrinos and antineutrino spectra obtained by solving the Boltzmann 
equation shows that, when the disk is in chemical equilibrium, the emission of neutrinos and 
antineutrinos are almost symmetric with nearly identical energy spectra, but the spectra is neithe r 
in the form of black body nor in the form of the gray body. And as shown in lPan &: Yuan! ([2012]), 
the black body spectra of neutrino a nd the neutrino spectr a based on the most commonly used 
simplified model of neutrino transport ()Di Matteo et al.ll2002l ) or equivalently the gray body spectra 
(jJaniuk et al.ll2007l ) can overestimate the annihilation luminosity by nearly one order of magnitude. 
In the following, we will also check the validity of the previous assumption on the neutrino transport. 

It is easy to estimate the mean optical depth of neutrinos in the inner part of the disk. From 
Fig. 6a and Fig. 2b, at r = 1 0M , the surface density is a = 3 x 10 17 g/cm 2 and the temperature 
is T = 7.4MeV. According to iDi Matteo et al.1 (|2002l ). the neutrino opacity of absorption r a and 
scattering r s are expressed as 

t s = 2.7T 2 1 cr 17 , (45) 

r a = 4.52^17, (46) 

where Tn is the temperature in unit of 10 11 K and oyj is the surface density in unit of 10 17 g/cm 2 . 
Therefore, the optical depth of neutrino s we get is t„ . = 10 , r s = 6, so the disk is optically thick 



for neutrinos at r = 10M. However, as iPan fc Yuan! (120121 ) has shown in quasi-o ptically opaque 
case (r as = 0.1 ~ 1), the neutrino spec tra are neither the black body spectra fhia ck (IPopham et al 



1993 ) nor the gray body spectra /gray (|Di Matteo et al 



2002 



■Taniuk et all 120071 ): 



/ black 



exp (p/kT) + 1 ' 



(47) 



Table 1: Luminosity of neutrino annihilation with different accretion rates 



M(M sun /s) 


T r=10M (MeV) 


L v + L p (ergs/s) 


L uP (ergs/s) 


V = L uP /(L u 


+ L D ) 


10. 


7.4 


1.04 x 10 54 


5.0 x 10 51 


4.8 x 10" 


-3 


1.0 


4.2 


0.92 x 10 53 


3.3 x 10 49 


3.6 x 10" 


-4 


0.1 


3.3 


0.30 x 10 52 


4.0 x 10 46 


1.3 x 10" 


-6 
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/gray 



exp (p/kT) + 1 ' 



(48) 



and the block factor b is given by 

b = 

where r = r a + r s , so here b = 0.16. 



(3/4)(r/2 + l/x/3 + l/3r a )' 



(49) 



For comparison, we plot the direction-averaged spectra F num of the black body spectra, the gray 
body spectra and the more realistic spectra obtained by solving Boltzmann equation (see Fig. 9). 
According to Fig. 9, it is obvious that the first assumption on neutrino transport overestimates the 
neutrino luminosity about 30% and the corresponding annihilation luminosity about (1.3 2 — 1) ~ 
70%, while the second assumption underestimates the neutrino luminosit y about 5 times and the 
corresponding annihilation luminosity about 25 times. It may explain whv lDi Matteo et al.l ((2002) 
drew a conclusion that the annihilation luminosity of NDAFs is no more than 1 50 e rgs/s, so 
NDAFs cannot serve as the central engine of GRBs, on the contrary. iPopham et aJj (|1999l ) claimed 
an opposite conclusion. 

The thickness of the disk is the third factor that affects the final annihilation luminosity: a 
larger ratio of thickness to radius means a larger solid angle the disk subtends to the annihilation 
zone. In the case we consider, the thick disk (H ~ r) will enhance the annihilation luminosity 
about 3 times than that of a thin disk. 

The signific ance of the distribution of electron fraction in NDAFs has been discussed by m any 
previous works (jKohri & Mineshige! \20od : \Lee et al.l bood ; Liu et al.l 120071 ; kohri et alJbood ). In 
this work, the distribution of electron fraction is a natural result of chemical equilibrium, thermal 
balance and hydrostatic equilibrium, rather than an artificial assumption in most previous works, 
so which should be most reliable. 




Fig. 9. — Comparison of the direction averaged black body spectra, gray body spectra and the real 
spectra obtained by solving the Boltzmann equation. 
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As shown in Fig. 5, helium is completely absent within 100M of NDAFs, so it is not needed 
to consider the contribution of helium to the total pressure and internal energy or the contribution 
of helium disintegration to the cooling of the disk. 

It should be emphasized that our calculation also has its own limitations. First, in this work, 
the disk is thick and has two dimensional structure, while the neutrino transport is treated in one 
dimension. Second, the effects of the motion of accretion flow and the curved spacetime on the 
vertical transport of neutrinos are neglected. These effects are significant in the inner part of the 
disk, for example, at r = 10M, the special relativity correction to the neutrino energy can be 
v/c ~ 30%, and the general relativity correction to the energy can be M/r ~ 10%. Third, in the 
zone near to the event horizon of the central black hole where the annihilation concentrates, the 
trajectories of neutrinos are severely bent by the central black hole and a considerable amount of 
neutrinos will undoubtedl y be captured by the hole. In addition, the gravitational instability of 



the disk was proposed by (jPerna et al.ll2006l ) to explain the energetic X-ray flares after the prompt 



emission in GRBs. The outer region of NDAFs with extremely high accretion rate is gravitational 



unstable (jChen &: Beloborodovll2007l ). will fragment and cause a variable accretion rate in the inner 
region. 
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